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Abstract 


A new three-dimensional mesh deformation algorithm, based on quaternion algebra, is 
introduced, A brief overview of quaternion algebra is provided, along with some preliminary 
results for two-dimensional structured and unstructured viscous mesh deformation. 


Introduction 

Mesh deformation is an important element in the analysis of moving bodies and shape 
optimization. The lack of robust and efficient mesh deformation tools is still a major barrier to 
routine applications of high-fidelity tools such as computational fluid dynamics (CFD) and 
computational structural mechanics (CSM) for multidisciplinary analysis and optimization. For 
example, CFD application for shape optimization requires a robust, automatic, and efficient tool 
to propagate the boundary deformation into the field mesh. For the gradient-based optimization, 
the efficiency is particularly crucial where in addition to the boundary deformation the sensitivity 
of the boundary coordinates must be propagated into the field mesh. Figure 1 shows an example 
of a boundary perturbation, where the boundary has been deformed, rotated, and translated. 
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Fig. 1 Undeformed and deformed Meshes 


The boundary deformation is defined as 

( 1 ) 

where r ; " are the undeformed boundary coordinates, and rf are the deformed boundary 
coordinates. There are two basic techniques to propagate the boundary perturbations into the 
field mesh: 1 ) mesh regeneration, and 2) mesh deformation. The next two subsections provide an 
overview of these techniques for structured and unstructured meshes. 


Structured Mesh 

Most structured grid regeneration and deformation techniques are based on transfinite 
interpolation (TFI). Gaitonde and Fiddes have provided a mesh regenerating technique based on 
TFI with exponential blending functions [l]. The choice ofblending functions has a considerable 
influence on the quality and robustness of the field mesh. Soni has proposed a set ofblending 
functions based on arclength [2]; such a set is extremely effective and robust for mesh 
regeneration and deformation. Jones and Samarehhave presented an algorithm for general 
multiblock mesh regeneration and deformation based on Soni's blending functions [3], 
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Hartwich and Agrawal have used a variation of the TFI method [4], They have introduced two 
new techniques: the use of the i slave-masteri concept to semiautomate the process, and the use 
of a Gaussian distribution function to preserve the integrity of meshes in the presence of multiple 
body surfaces. Wong et al. have used Algebraic and Iterative Mesh 3D (AIM3D), which is based 
on a combination of algebraic and iterative methods [5], Leatham and Chappell have used a 
Laplacian technique more commonly used for unstructured mesh deformation [6], 

Unstructured Mesh 

For unstructured meshes with large geometry changes, a new mesh may need to be regenerated 
at the beginning of each optimization cycle. Botkin has introduced a local remeshing procedure 
that operates only on the specific edges and faces associated with the design variable changes [7], 
Similarly, Kodiyalam, Kumar, and Finnigan have used a mesh regeneration technique based on 
the assumption that the solid model topology stays fixed for small perturbations [8], Solid model 
topology comprises the number of mesh-points, edges, and faces. Any change in the topology 
will cause the model regeneration to fail. To avoid such a failure, a set of constraints among 
design variables must be satisfied, in addition to constraints on their bounds. 

During shape optimization, the boundary mesh may undergo many small deformations; it would 
be too costly to regenerate the mesh in response to these deformations. In addition, the new, 
regenerated mesh may not have the same number of mesh points and/or the same connectivity. 
Either of these situations will result in discontinuous sensitivity derivatives. Batina has presented 
a mesh deformation algorithm that did alleviate the need for mesh regeneration. Batinais 
approach models mesh edges with springs [9]. The spring stiffness k, k for a given edge jk is 

taken to be inversely proportional to the element edge length. Then, the field mesh movement is 
computed through the static equilibrium equations: 
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where k Jk = 



( 2 ) 


The summation is over all the edges of the elements. The coefficient k jk is relatively large for 

small cells. Therefore these small cells, which are usually near the boundary of the body, tend to 
undergo rigid body movement. This rigid body movement avoids rapid variations in 
deformation, thus eliminating the possibility of small cells having very large changes in volume. 
These large changes could lead to negative cell volumes. 


Blom [10] has provided a detailed analysis for the spring method and draws an analogy between 
the spring method and an elliptic differential equation approach for structured mesh generation. 
Zhang and Belegundu have proposed an algorithm similar to the spring analogy that can handle 
large mesh deformation [ll]. They have used the ratio of the cell Jacobian to the cell volume for 
the spring stiffness. Crumpton and Giles have found the spring analogy inadequate and 
ineffective for large mesh deformations [ 12] and proposed a formulation based on the heat 
conduction equation with the coefficient of thermal conductivity inversely proportional to cell 
volume. They attributed their success to the choice of cell volume used in the criteria for a valid 
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mesh. In contrast, the spring analogy uses only edges, which are not directly linked to the mesh 
validity. 

Farhat et al. [13] have proposed a modification to the spring analogy algorithm to include 
additional torsional spring to control mesh skewness and folding. For two-dimensional 
applications, they demonstrated that the modified algorithm has advantages in terms of 
robustness, quality, and performance. 

Tezduyar and Behr [14] have proposed an algorithm based on linear elasticity, which includes 
full stress tensor. Cavallo et al. [15] have applied this method to mesh deformation for 
aero/propulsive flowfield calculations. They noted that the method preserves the mesh quality, 
and it produces a better mesh than the spring analogy method. The linear elasticity approach 
requires solving the complete stress tensor. In contract, the spring analogy represents only the 
diagonal elements of the stress tensor. Cavallo et al. have concluded that the elasticity approach 
is considerably more expensive. 


Role of boundary orientation in mesh deformation 

The traditional deformation algorithms, such as interpolation and spring analogy, use boundary 
translation to deform the field mesh. However, the boundary deformation alters the boundary 
position as well as the boundary orientation (i.e., rotation angle) as shown in Fig. (2). The 
traditional mesh deformation algorithms do not use this additional information on the changes in 
the boundary orientation. Morton, Melville, and Visbal [16] have proposed a TFI algorithm to 
interpolate the boundary deformation as well as the changes in the orientation through Euler 
angle. They concluded the inclusion of Euler angle preserves the mesh orthogonality for 
significant deformations. They successfully applied the algorithm to a two-dimensional 
structured CFD mesh. 



Extension of the approach of Morton, Melville, and Visbal [16] to three-dimensional applications 
requires the direct interpolation of the changes in the orientation (three Euler angles). 

Historically, Euler angle representation is the most popular interpolation technique for 
orientation [17], Euler angles ignore the interaction of multiple rotations about the separate axes 
which could lead to i gimbal lock! as described by Watt and Watt [17], 
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The three Euler rotations can be accomplished by a single rotation about a vector. This single 
rotation simplifies the interpolation process, but it has the inherent problem of non-smooth 
interpolation and the so-called i gimbal-lock.i To avoid these problems, quaternions are used to 
represent the changes in the boundary orientation. 

A brief introduction of quaternion algebra is presented in the next section. Then, a general three- 
dimensional mesh deformation algorithm based on quaternion algebra is presented. 


What are quaternions? 


Only brief review of quaternion algebra is provided here; readers are referred to the work of 
Altmann [18], Shoemake [19], and Philips, Hailey, and Gerbert [20] for more details. There is 
some controversy on who invented quaternion algebra. The articles by Altman [18] and Philips 
et al. [20] provide a very interesting history of quaternion algebra. 

A quaternion is a generalized complex number (hypercomplex number) that is composed of one 
real and three imaginary numbers (0 = q 0 + q x i +q 2 j + q i k ) , where ii = jj = kk = - 1, 
ij = -ji = k, jk = -kj = i, ki = -ik = j . The following is a set quaternion properties that will be 
used later: 

• Conjugate of a quaternion, Q* = q 0 - qq - q 2 j -q 3 k 

• Magnitude of a quaternion, ||^|| = ^00 = ■yjql +qf +ql +q] 

• Unit quaternion, ||0|| = 1 

• Associative, (Q ] Q 2 )Q 3 = Q X {Q 2 Q 3 ) 

• Not commutative, Q X Q 2 * 0 1 0 [ 

• Inverse of a quaternion, Q~ x = Q* KQQ* ) 

• For unit quaternions, Q~' =Q * 

A quaternion can be inteipreted as a scalar together with a vector (direction), 

Q = [s,\],s = q 0 ,\ = (q\,q 2 ,q 3 ) 

In this notation, quaternion multiplication has the particularly simple form 
Q X Q 2 =[^,v 1 ][^ 2 ,v 2 ] = [^ 1 5 2 - v, • v 2 ,JjV 2 +s 2 v, + v, xv 2 ] 

where • denotes the vector dot product, and x denotes the vector cross product. 

Quaternions are ideal for modeling rotations. The last three components of a quaternion represent 
the axis around which the rotation occurs, and the first component represents the magnitude of 
the rotation. There are three steps involved in rotating a point, p , about a unit vector, u , by an 
angle, 0 . First, a quaternion is constructed for the point as P = [0,p] . Second, a quaternion is 
constructed for the rotation as 


Q 


r i 0 • e 

= U, vj,^ = cos— , v = usin — 
2 2 


( 3 ) 


Third, the point is rotated as P mlillcA = QPQ 1 . If Q is a unit quaternion, then we can use the 


conjugate of the quaternion to perform the rotation, P mtated = QPO* . Multiple rotations can be 
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simplified by using a single quaternion. For example, if Q l and Q 2 are unit quaternions 
representing two rotations, the two rotations can be combined as 

Q 2 (Qm [ )Q~2 = (Q 2 q ] )P( 0; ] 0 2 ') = (Q 2 Q)P(Q 2 Q,y' 

Quaternion coordinates represent rotation as Cartesian coordinates represent translation as a 
single vector. This characteristic has been fully exploited in representing attitude of aircraft 
kinematics [20], Quaternion coordinates are best for interpolation of orientation as used in 
computer animation. Shoemake has presented a robust and efficient application of quaternions 
for BEzier interpolation of orientation used in computer animation [19]. 


Quaternions and Mesh Deformation 

This section presents a technique to model the boundary deformation by quaternion algebra. 
When these boundary quaternions are applied to the undeformed boundary mesh, they produce 
the deformed boundary mesh and orientation. 

The deformation vectors, 8 , represent the boundary translation, which is defined in the Euclidian 
space. In traditional mesh deformation algorithms, these vectors are used to propagate the 
deformation into the field mesh. In a similar manner, we will use the boundary quaternions to 
propagate the deformation into the field mesh. The process of determining boundary quaternions 
is divided into three steps, as shown in Fig. (3). 



Fig. 3 Process of boundary quaternion construction 
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In step 1 , the mesh points for the undeformed boundary ( r“ ), the deformed boundary ( rf ), and 
the neighboring points ( r" and r d ) are translated to the origin. 


U d 1 d 

■ r t , r = r 


(4) 


In the second step, r Ml is rotated so that the undeformed boundary normal vector aligns with the 
deformed boundary normal vector. This rotation is modeled with a quaternion. First, the normal 
vector of a plane shared by both deformed and undeformed normal vectors share ( defined as 
u = n a x n ) and the angle a between two normal vectors is determined. Then, a quaternion is 


defined for the rotation as Q x = [cos ^ > n x n</ sin . Points r" 1 are rotated by quaternion to 


form r" 2 , such that 


r a2 =Q[O,r" 1 ]0 1 - 1 (5) 

In the third step, points r“ 2 are rotated about the deformed boundary normal vector to minimize 
the angle between corresponding neighboring points. The optimum rotation angle, 6 , is defined 
as the average angle between corresponding edges of r a2 and the edges of deformed boundary. 
Another quaternion can then be defined for this rotation, Q 2 = [cos ^ - n<l s ' n %] • 


These two quaternions are combined to form a single quaternion as Q t = Q^Q 2 . The total 
translation vector for the boundary can now be defined as A, = rf - d,., where 

[d,.,0] = a[0,r i d ]a.- 1 . 

Quaternions and total translation vectors for all boundary mesh points have been computed. The 
translation vectors account for the translation, and quaternions account for the changes in the 
boundary orientation. 

The translation vectors and quaternions are propagated into the field mesh by one of the 
traditional deformation algorithms such as TFI or the spring analogy. Then, the field mesh is 
updated based on the field values for the translation vectors and quaternions as 

Afield = Afield + QfieldW’^field^Q field (Q 

Results 

The results are presented for structured and unstructured viscous mesh deformations. Figure (4) 
shows a viscous structured mesh with 257x65 mesh points. The undeformed mesh lines are 
orthogonal to the boundaries. The boundary mesh is deformed, rotated, and translated to simulate 
aeroelastic deformation. Figure (4) shows comparisons of TFI (left side of figure) and quaternion 
approach (right side of figure). Unlike the traditional TFI, the quaternion approach can clearly 
preserve the boundary orthogonality. Because the boundary quaternions are based on the changes 
in the boundary mesh point positions as well as the orientations, the algorithm can guarantee that 
the mesh near the boundary has the same characteristics as the undeformed mesh. Figure (4) 
clearly demonstrates this important property. 
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Next, the quaternion approach is applied to an unstructured viscous mesh, where the flap has 
been rotated. The spring analogy was used to propagate the boundary quaternions to the field 
mesh. The results are shown in Fig. (5). Again, the use of quaternion has preserved the mesh 
characteristics. 

Conclusions 

A new three-dimensional mesh deformation algorithm based on quaternion algebra has been 
presented. These preliminary two-dimensional results indicate the traditional algorithms such as 
TFI and spring analogy can be easily augmented with the quaternions to preserve mesh quality 
near the viscous boundary. We plan to apply this method for three-dimensional structured and 
unstructured viscous meshes. We also plan to evaluate the quality of meshes deformed by 
quaternion approach by means of CFD applications. 
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Fig. 4 Deformation comparison for structured mesh 



Fig. 5 Deformation comparison for unstructured mesh 
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